# SOURCE: https://cedricscherer.netlify.app/2019/05/17/the-evolution-of-a-ggplot-ep.-1/#aim 

# Packages
required_packages <- c(
                       "tidyverse",   
                       "rgdal",
                       "RColorBrewer",
                       "leaflet",
                       "rgeos",
                       "openxlsx",
                       "ggthemes",
                       "tidyverse"
                       )      

for (pkg in required_packages) {
  # install packages if not already present
  if (!pkg %in% rownames(installed.packages())) {
    install.packages(pkg)
  }
  
  # load packages to this current session 
  library(pkg, character.only = TRUE)
}

library(tidyverse)

df_students <- readr::read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2019/2019-05-07/student_teacher_ratio.csv")

df_world_tile <- readr::read_csv("https://gist.githubusercontent.com/maartenzam/787498bbc07ae06b637447dbd430ea0a/raw/9a9dafafb44d8990f85243a9c7ca349acd3a0d07/worldtilegrid.csv") %>% 
  mutate(
    ## Namibias two-digit country code is handled as `NA` - let us fix that
    alpha.2 = if_else(name == "Namibia", "NA", alpha.2),
    ## We are going to split "Americas" into "North America" and "Sout America"
    region = if_else(region == "Americas", sub.region, region),
    region = if_else(region %in% c("Northern America", "Central America", "Caribbean"), 
                     "North America", region),
    region = if_else(region == "Southern America", "South America", region),
    ## to join both data sets, we need a id column
    country_code = alpha.3
  )

df_ratios <- df_students %>% 
  ## Let's keep only the most recent data per country
  group_by(country, indicator) %>% 
  filter(year == max(year)) %>% 
  ungroup() %>% 
  # Create `NA`s for countries which do not have any data 2012-2018
  complete(indicator, nesting(country, country_code)) %>% 
  ## Let's focus on primary education and keep only countries (coded by letters)
  filter(
    indicator == "Primary Education",
    str_detect(country_code, "[A-Z]")
  ) %>% 
  ## merge with world tile map data
  full_join(df_world_tile) %>%
  filter(
    !is.na(region),
    !is.na(indicator)
  ) %>% 
  group_by(region) %>% 
  mutate(student_ratio_region = median(student_ratio, na.rm = T)) %>% 
  ungroup()

p1<- ggplot(df_ratios, aes(x = region, y = student_ratio)) +
  geom_boxplot()
p1

# Sort the data
df_sorted <-
  df_ratios %>%
  mutate(region = fct_reorder(region, -student_ratio_region))


p2<- ggplot(df_sorted, aes(x = region, y = student_ratio)) +
  geom_boxplot()
p2

p2+ labs(title = "Student to teacher ratio",
         subtitle = "Lastest available data",
         y = "Student to teacher ratio",
         x = "Region", 
         caption = "Source: UNESCO")

p3<- ggplot(df_sorted, aes(x = region, y = student_ratio)) +
  geom_boxplot() +
  coord_flip() +
  scale_y_continuous(limits = c(0, 90)) +
  labs(title = "Student to teacher ratio",
         subtitle = "Lastest available data",
         y = "Student to teacher ratio",
         x = "Region", 
         caption = "Source: UNESCO")

p3

ggsave( filename = here::here("figures", 
                  paste0(Sys.Date(), "_hospital_debt", ".png")),
                  width = 20.49,
                  height = 10.3,
                  units = c("cm"),
                  dpi = "retina")

comm_pt <- readOGR(
  dsn= "datasets/concelhos-shapefile" , ###
  layer="concelhos",
  verbose=FALSE,
  use_iconv = TRUE,   ###  
  encoding = "UTF-8"  ### 
)

# make the polygons a bit less verbose
comm_pt1 <- gSimplify(comm_pt, 0.01, topologyPreserve=TRUE)
comm_pt1 = SpatialPolygonsDataFrame(comm_pt1, data=comm_pt@data)

comm_pt1@data <- comm_pt1@data %>%
  dplyr::select(
    "ISO",
    "ID_1",
    "NAME_1",
    "ID_2",
    "NAME_2",
  )

covid_inc <- read_csv("datasets/covid_inc.csv") 


covid_inc <- covid_inc %>%
  rename(NAME_2=Concelho) %>%
  mutate( NAME_2=recode(NAME_2,
                         'Ponte de Sor' = "Ponte de Sôr"
                         ))


comm_pt1@data<-left_join(comm_pt1@data,
                     covid_inc,
                     by = "NAME_2")



#comm_pt1@data<- comm_pt1@data %>%
#  filter(
#    !ARS %in% c("Açores", "Madeira")
#  )

mybins <- c(0,60,120,240,480,Inf)
mycolors <- c("#f7f7f7","#fee391","#fb6a4a","#cb181d","#67000d")
mypalette <- colorBin( palette=mycolors, domain=comm_pt1@data$Incidencia, na.color="transparent", bins=mybins)

# Prepare the text for tooltips:
mytext <- paste(
    "Concelho: ", comm_pt1@data$NAME_2,"<br/>",
    "Incidência: ", comm_pt1@data$Incidencia, "<br/>",
    "Categoria: ", comm_pt1@data$Incidencia_descritivo,
    sep="") %>%
  lapply(htmltools::HTML)

# Final Map
leaflet(comm_pt1) %>%
  addTiles()  %>%
  setView( lat=39.5, lng=-7.6, zoom=7) %>%
  addProviderTiles("CartoDB.Positron") %>%
  addPolygons(
    fillColor = ~mypalette(Incidencia),
    stroke=TRUE,
    fillOpacity = 0.9,
    color="grey",
    weight=0.3,
    label = mytext,
    labelOptions = labelOptions(
      style = list("font-weight" = "normal", padding = "3px 8px"),
      textsize = "13px",
      direction = "auto"
    )
  ) %>%
  addLegend( pal=mypalette, values=~Incidencia, opacity=0.9, title = "Incidência", position = "bottomleft" )

library(viridis)

# I need to fortify the data AND keep trace of the commune code! (Takes ~2 minutes)
library(broom)
comm_pt1_fortified <- tidy(comm_pt1, region = "NAME_2")

ggplot() +
  geom_polygon(data=comm_pt1_fortified,
               aes(x=long,
                   y=lat,),
               color="white",
               size=.2) +
  theme_void() +
  coord_map() +
  scale_fill_binned(type = "viridis", direction=-1) +
  labs(title="Population in New York City",
       subtitle="Neighborhoods are filled by population",
       fill="Population")